#install.packages("pacman")#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)#devtools::install_github("daniel1noble/orchaRd", force = TRUE)library(tidyverse)library(here)library(lme4)library(orchaRd)#library(gptstudio)library(metafor)library(patchwork)library(alluvial)library(ggalluvial)library(easyalluvial)library(ape)library(clubSandwich)library(emmeans)library(MuMIn)# making metafor talk to MuMIneval(metafor:::.MuMIn)# install.packages("pak")#pak::pak("MichelNivard/gptstudio")
2 Getting data loaded
Code
#dat_full <- read.csv(here("data/dat_04_04_2023.csv"))dat_full <-read.csv(here("data/dat_28_06_2023.csv"))# add phylogenetic tree - only topologies# TODO? - we could get better tree from birdtree.org# we can do 50 different trees as in # https://academic.oup.com/sysbio/article/68/4/632/5267840tree_top <-read.tree(here("R/birds_MA.tre"))# tree with branch lengthstree <-compute.brlen(tree_top)#plot(tree)# turning it into a correlation matrixcor_tree <-vcv(tree,corr=T)
3 Preparing data set
Code
# function for calculating varianceVd_func <-function(d, n1, n2, design, r =0.5){# independent designif(design =="among"){ var <- (n1 + n2) / (n1*n2) + d^2/ (2* (n1 + n2 -2)) # variance } else { # dependent design var <-2*(1-r) / n1 + d^2/ (2*(n1 -1)) # variance } var # return variance}# getting Hedges' g - get small size bias corrected effect sizedat_full$SMD <- dat_full$d / (1-3/(4* (dat_full$NTreat + dat_full$Ncontrol) -9))# flipping d dat_full$SMD <- dat_full$d*dat_full$Direction*dat_full$PredictedDirection# calucating Vddat_full$Vd <-with(dat_full, pmap_dbl(list(SMD, NTreat, Ncontrol, Design), Vd_func))# extra useful function# function for getting mean and sd from median, quartiles and sample size# get_mean_sd <- function(median, q1, q3, n){# sd <- (q3 - q1) / (2 * (qnorm((0.75 * n - 0.125) / (n + 0.25)))) # sd# mean <- (median + q1 + q3)/3 # mean# c(mean, sd)# }# observation iddat_full$Obs_ID <-1:nrow(dat_full)dat_full$Phylo <-gsub(" ", "_", dat_full$FocalSpL)# filtering very large variance and also very small sample sizedat_int <- dat_full %>%filter(Vd <10& Ncontrol >2& NTreat >2)#dim(dat_full)#dim(dat_int)# sorting out modality stuff# creat - 1,2,3 modality - also easier classification A, O, V (AOV = L) dat_int %>%mutate(Treat_mod =case_when(Treatment =="A"~"A", Treatment =="AV"~"AV", Treatment =="AVG"~"AV", Treatment =="AVM"~"AV", Treatment =="L"~"AVO", Treatment =="O"~"O", Treatment =="OV"~"OV", Treatment =="V"~"V", Treatment =="VG"~"V", Treatment =="VM"~"V", Treatment =="VP"~"V"),# into how manyTreat_No =case_when(Treatment =="A"~1, Treatment =="AV"~2, Treatment =="AVG"~2, Treatment =="AVM"~2, Treatment =="L"~3, Treatment =="O"~1, Treatment =="OV"~2, Treatment =="V"~1, Treatment =="VG"~1, Treatment =="VM"~1, Treatment =="VP"~1),# des it have some add-onsAdd_on =case_when(Treatment =="A"~"No", Treatment =="AV"~"No", Treatment =="AVG"~"Yes", Treatment =="AVM"~"Yes", Treatment =="L"~"No", Treatment =="O"~"No", Treatment =="OV"~"No", Treatment =="V"~"No", Treatment =="VG"~"Yes", Treatment =="VM"~"Yes", Treatment =="VP"~"Yes"), ) -> dat# creating data just for A, V, and AV dat_short <- dat %>%filter(Treat_mod =="A"| Treat_mod =="V"| Treat_mod =="AV")# for add-on, we only need V and AVdat_short_add <- dat %>%filter(Treat_mod =="AV"| Treat_mod =="V")dat <- dat %>%mutate_if(is.character, as.factor)#summary(dat)
4 Exploratory visualization
Code
# some data exploration# dat %>% group_by(Treat_mod) %>% summarise(n = n())# using a table to see overalps between Treat_mod and Typetab <-table(dat_short$Treat_mod, dat_short$Type)# visualise this table using alluvial plot?# https://cran.r-project.org/web/packages/alluvial/vignettes/alluvial.html# Treatment vs Typedat_short %>%group_by(Treat_mod, Type) %>%summarise(n =n()) -> tab1#alluvial(tab1[,1:2], freq = tab1$n)# using ggaruvialggplot(tab1,aes(y = n,axis1 = Treat_mod,axis2 = Type)) +geom_alluvium(aes(fill = Treat_mod)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =6, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Treatment modality and trait type")
Code
# Type vs Sexdat %>%group_by(Sex, Type) %>%summarise(n =n()) -> tab2#alluvial(tab2[,1:2], freq = tab2$n)# using ggaruvialggplot(tab2,aes(y = n,axis1 = Type,axis2 = Sex)) +geom_alluvium(aes(fill = Type)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =6, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Trait type and sex")
Code
# TOOD - use easyalluvial to make one fig# other ones - easyalluvial#https://www.r-bloggers.com/2018/10/data-exploration-with-alluvial-plots-an-introduction-to-easyalluvial/# using easyalluvial and alluvial_wide#factor_cols <- dat %>% select_if(is.factor) %>% names()#alluvial_wide(dat_short, , max_variables = 5# , fill_by = 'first_variable' ) %>%# add_marginal_histograms(dat_short)# exploratory analysis# check each columns for missing values and other stuff# dat %>% map_df(~sum(is.na(.)))
orchard_plot(mod1, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
6.2 Treatment with relevant data
Code
VCVs <-vcalc(vi = dat_short$Vd,cluster = dat_short$SubjectID,rho =0.5)mod1a <-rma.mv(yi = SMD,V = VCVs,random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Treat_mod, test ="t",method ="REML", sparse =TRUE,data = dat_short)# to look at A vs AV and A vs V# TODO - make these hetero model toosummary(mod1a)
mod1b <-rma.mv(yi = SMD, V = VCVs, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~relevel(factor(Treat_mod), ref ="AV"), test ="t",method ="REML", sparse =TRUE,data = dat_short)# to look at AV vs V and AV vs A (already done)summary(mod1b)
df AIC BIC AICc logLik LRT pval QE
Full 8 1663.1371 1698.1237 1663.3867 -823.5686 NA
Reduced 6 1680.3773 1706.6172 1680.5224 -834.1887 21.2402 <.0001 NA
Code
orchard_plot(mod1c, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
6.3 Treatment with some additions
Code
# the effect of additions# this is a part of sensitivity analysisVCVs2 <-vcalc(vi = dat_short_add$Vd,cluster = dat_short_add$SubjectID,rho =0.5)mod5 <-rma.mv(yi = SMD, V = VCVs2, random =list(~1|FocalSpL , ~1| RecNo, ~ Treat_mod | Obs_ID), mod =~ Treat_mod*Add_on -1,test ="t",struct ="DIAG",method ="REML", sparse =TRUE,data = dat_short_add)summary(mod5)
# testing the number of stimulimod4 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Treat_No, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod4)
bubble_plot(mod4,mod ="Treat_No",group ="RecNo",xlab ="The number of simuli",g =TRUE)
6.5 Trait type
Code
# Type of responsesmod2 <-rma.mv(yi = SMD, V = VCV, random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Type -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod2)
# heteroscadasticity model better than the homoscedasticity model# note LifeHistory has lowest variation but this may be expected? # as it is less flexiable (e.g. the number of eggs?)anova(mod2, mod2b)
df AIC BIC AICc logLik LRT pval QE
Full 8 1722.2085 1757.7233 1722.4419 -853.1043 NA
Reduced 6 1757.0783 1783.7144 1757.2140 -872.5392 38.8698 <.0001 NA
Code
orchard_plot(mod2b, mod ="Type",group ="RecNo",xlab ="Standardised mean differnece (SMD)",branch.size =3)
6.6 trait categories
Code
# Category of responsesmod3 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Category -1, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod3)
orchard_plot(mod3, mod ="Category",group ="RecNo", xlab ="Standardised mean differnece (SMD)",angle =45,branch.size =3)
6.7 Predactor guild
Code
# Predactor guild# quite heterogeneous# TODO this could be in random effects - think abou thtis a bit latermod6 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ PredGuild -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod6)
orchard_plot(mod11,mod ="ControlType",group ="RecNo",xlab ="Standardised mean differnece (SMD)")
6.13 Sex
Code
# sex# TODO - this could be interesting# what is in males and femalesmod12 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Sex -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod12)
bubble_plot(mod15,mod ="ln_duration",group ="RecNo",xlab ="log(duration in days)",g =TRUE) +geom_point(data = dat,aes(x = ln_duration, y = SMD,color = Type,fill = Type,size =1/sqrt(Vd)), alpha =0.6) +scale_color_discrete() +#+ # how to put the legend for colourguides(color ="legend")
orchard_plot(mod_full,mod ="Type",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
orchard_plot(mod_full,mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
int_type <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Type")bubble_plot(int_type, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3)
Code
int_trt <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Treat_mod")bubble_plot(int_trt, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3)
# displays delta AICc <2candidates_aic2 <-subset(candidates, delta <5) # model averagingmr_averaged_aic2 <-summary(model.avg(candidates, delta <5)) # relative importance of each predictor for all the modelsimportance <-sw(candidates)
dat_fulle <-qdrg(object = mod_fulle, data = dat_short)# marginalized overall mean at vi = 0 and year.c = 0; also weights = "prop" or "cells" average things over proportionally. if not specified, all groups (levels) get the same weights# res_fulle1 <- emmeans(dat_fulle, # specs = ~ sqeffectN,# df = mod_fulle$ddf, # weights = "prop")
Source Code
---title: "**Integration of multimodal cues does not alter mean but reduces variance in avian responses to predators: a systematic review and meta-analysis**"author: "**Kim + Shinichi et al.**"date: "`r Sys.Date()`"format: html: toc: true toc-location: left toc-depth: 3 toc-title: "**Table of Contents**" output-file: "index.html" theme: simplex embed-resources: true code-fold: true code-tools: true number-sections: true #bibliography: ./bib/ref.bib fontsize: "12" max-width: "10" code-overflow: wrapcrossref: fig-title: Figure # (default is "Figure") tbl-title: Table # (default is "Table") title-delim: — # (default is ":") fig-prefix: Fig. # (default is "Figure") tbl-prefix: Tab. # (default is "Table")editor_options: chunk_output_type: consoleexecute: warning: false message: false tidy: true cache: true---## Setting up```{r}#install.packages("pacman")#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)#devtools::install_github("daniel1noble/orchaRd", force = TRUE)library(tidyverse)library(here)library(lme4)library(orchaRd)#library(gptstudio)library(metafor)library(patchwork)library(alluvial)library(ggalluvial)library(easyalluvial)library(ape)library(clubSandwich)library(emmeans)library(MuMIn)# making metafor talk to MuMIneval(metafor:::.MuMIn)# install.packages("pak")#pak::pak("MichelNivard/gptstudio")```## Getting data loaded```{r}#dat_full <- read.csv(here("data/dat_04_04_2023.csv"))dat_full <-read.csv(here("data/dat_28_06_2023.csv"))# add phylogenetic tree - only topologies# TODO? - we could get better tree from birdtree.org# we can do 50 different trees as in # https://academic.oup.com/sysbio/article/68/4/632/5267840tree_top <-read.tree(here("R/birds_MA.tre"))# tree with branch lengthstree <-compute.brlen(tree_top)#plot(tree)# turning it into a correlation matrixcor_tree <-vcv(tree,corr=T)```## Preparing data set```{r}# function for calculating varianceVd_func <-function(d, n1, n2, design, r =0.5){# independent designif(design =="among"){ var <- (n1 + n2) / (n1*n2) + d^2/ (2* (n1 + n2 -2)) # variance } else { # dependent design var <-2*(1-r) / n1 + d^2/ (2*(n1 -1)) # variance } var # return variance}# getting Hedges' g - get small size bias corrected effect sizedat_full$SMD <- dat_full$d / (1-3/(4* (dat_full$NTreat + dat_full$Ncontrol) -9))# flipping d dat_full$SMD <- dat_full$d*dat_full$Direction*dat_full$PredictedDirection# calucating Vddat_full$Vd <-with(dat_full, pmap_dbl(list(SMD, NTreat, Ncontrol, Design), Vd_func))# extra useful function# function for getting mean and sd from median, quartiles and sample size# get_mean_sd <- function(median, q1, q3, n){# sd <- (q3 - q1) / (2 * (qnorm((0.75 * n - 0.125) / (n + 0.25)))) # sd# mean <- (median + q1 + q3)/3 # mean# c(mean, sd)# }# observation iddat_full$Obs_ID <-1:nrow(dat_full)dat_full$Phylo <-gsub(" ", "_", dat_full$FocalSpL)# filtering very large variance and also very small sample sizedat_int <- dat_full %>%filter(Vd <10& Ncontrol >2& NTreat >2)#dim(dat_full)#dim(dat_int)# sorting out modality stuff# creat - 1,2,3 modality - also easier classification A, O, V (AOV = L) dat_int %>%mutate(Treat_mod =case_when(Treatment =="A"~"A", Treatment =="AV"~"AV", Treatment =="AVG"~"AV", Treatment =="AVM"~"AV", Treatment =="L"~"AVO", Treatment =="O"~"O", Treatment =="OV"~"OV", Treatment =="V"~"V", Treatment =="VG"~"V", Treatment =="VM"~"V", Treatment =="VP"~"V"),# into how manyTreat_No =case_when(Treatment =="A"~1, Treatment =="AV"~2, Treatment =="AVG"~2, Treatment =="AVM"~2, Treatment =="L"~3, Treatment =="O"~1, Treatment =="OV"~2, Treatment =="V"~1, Treatment =="VG"~1, Treatment =="VM"~1, Treatment =="VP"~1),# des it have some add-onsAdd_on =case_when(Treatment =="A"~"No", Treatment =="AV"~"No", Treatment =="AVG"~"Yes", Treatment =="AVM"~"Yes", Treatment =="L"~"No", Treatment =="O"~"No", Treatment =="OV"~"No", Treatment =="V"~"No", Treatment =="VG"~"Yes", Treatment =="VM"~"Yes", Treatment =="VP"~"Yes"), ) -> dat# creating data just for A, V, and AV dat_short <- dat %>%filter(Treat_mod =="A"| Treat_mod =="V"| Treat_mod =="AV")# for add-on, we only need V and AVdat_short_add <- dat %>%filter(Treat_mod =="AV"| Treat_mod =="V")dat <- dat %>%mutate_if(is.character, as.factor)#summary(dat)```## Exploratory visualization```{r}# some data exploration# dat %>% group_by(Treat_mod) %>% summarise(n = n())# using a table to see overalps between Treat_mod and Typetab <-table(dat_short$Treat_mod, dat_short$Type)# visualise this table using alluvial plot?# https://cran.r-project.org/web/packages/alluvial/vignettes/alluvial.html# Treatment vs Typedat_short %>%group_by(Treat_mod, Type) %>%summarise(n =n()) -> tab1#alluvial(tab1[,1:2], freq = tab1$n)# using ggaruvialggplot(tab1,aes(y = n,axis1 = Treat_mod,axis2 = Type)) +geom_alluvium(aes(fill = Treat_mod)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =6, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Treatment modality and trait type")# Type vs Sexdat %>%group_by(Sex, Type) %>%summarise(n =n()) -> tab2#alluvial(tab2[,1:2], freq = tab2$n)# using ggaruvialggplot(tab2,aes(y = n,axis1 = Type,axis2 = Sex)) +geom_alluvium(aes(fill = Type)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =6, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Trait type and sex")# TOOD - use easyalluvial to make one fig# other ones - easyalluvial#https://www.r-bloggers.com/2018/10/data-exploration-with-alluvial-plots-an-introduction-to-easyalluvial/# using easyalluvial and alluvial_wide#factor_cols <- dat %>% select_if(is.factor) %>% names()#alluvial_wide(dat_short, , max_variables = 5# , fill_by = 'first_variable' ) %>%# add_marginal_histograms(dat_short)# exploratory analysis# check each columns for missing values and other stuff# dat %>% map_df(~sum(is.na(.)))```## Meta-analysis### All random effects```{r}##| warning: false# VCV matrixVCV <-vcalc(vi = dat$Vd,cluster = dat$SubjectID,rho =0.5)mod0 <-rma.mv(yi = SMD,V = VCV, random =list(~1| Phylo,~1| FocalSpL,~1| RecNo,~1| SubjectID, # incoprated as VCV~1| Obs_ID),R =list(Phylo = cor_tree),test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod0)# TODO - think about whether we add this or notrobust(mod0, cluster = dat$SubjectID)round(i2_ml(mod0), 2)orchard_plot(mod0,group ="RecNo",xlab ="Standardised mean differnece (SMD)",branch.size =3)```### Reduced model```{r}# reduced modelmod0r <-rma.mv(yi = SMD,V = VCV, random =list(#~1 | Phylo,~1| FocalSpL,~1| RecNo,#~1 | SubjectID, # incoprated as VCV~1| Obs_ID),#R = list(Phylo = cor_tree),test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod0r)round(i2_ml(mod0r), 2)```## Meta-regression: uni-moderator (mostly)### Treatmeant with all data```{r}## Treatment - A, V, AV etc #dat$Phylo <- as.character(dat$Phylo)#match(dat$Phylo, colnames(cor_tree))#match(colnames(cor_tree), dat$Phylo)mod1 <-rma.mv(yi = SMD, V = VCV, random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Treat_mod -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod1)round(r2_ml(mod1)*100, 2)orchard_plot(mod1, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)```### Treatment with relevant data```{r}VCVs <-vcalc(vi = dat_short$Vd,cluster = dat_short$SubjectID,rho =0.5)mod1a <-rma.mv(yi = SMD,V = VCVs,random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Treat_mod, test ="t",method ="REML", sparse =TRUE,data = dat_short)# to look at A vs AV and A vs V# TODO - make these hetero model toosummary(mod1a)mod1b <-rma.mv(yi = SMD, V = VCVs, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~relevel(factor(Treat_mod), ref ="AV"), test ="t",method ="REML", sparse =TRUE,data = dat_short)# to look at AV vs V and AV vs A (already done)summary(mod1b)orchard_plot(mod1a, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)# modeling heteroscedasticitymod1c <-rma.mv(yi = SMD, V = VCVs, random =list(~1|FocalSpL , ~1| RecNo, ~ Treat_mod | Obs_ID), mod =~ Treat_mod, test ="t",struct ="DIAG",method ="REML", sparse =TRUE,data = dat_short)summary(mod1c)# comparision modelsanova(mod1a, mod1c)orchard_plot(mod1c, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)```### Treatment with some additions```{r}# the effect of additions# this is a part of sensitivity analysisVCVs2 <-vcalc(vi = dat_short_add$Vd,cluster = dat_short_add$SubjectID,rho =0.5)mod5 <-rma.mv(yi = SMD, V = VCVs2, random =list(~1|FocalSpL , ~1| RecNo, ~ Treat_mod | Obs_ID), mod =~ Treat_mod*Add_on -1,test ="t",struct ="DIAG",method ="REML", sparse =TRUE,data = dat_short_add)summary(mod5)```### Treatment as an ordinal variable```{r}# testing the number of stimulimod4 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Treat_No, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod4)bubble_plot(mod4,mod ="Treat_No",group ="RecNo",xlab ="The number of simuli",g =TRUE)```### Trait type```{r}# Type of responsesmod2 <-rma.mv(yi = SMD, V = VCV, random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Type -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod2)mod2c <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Type,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod2c)mod2d <-rma.mv(yi = SMD,V = VCV,random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mod =~relevel(Type, ref ="LifeHistory"),test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod2d)orchard_plot(mod2,mod ="Type",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)round(r2_ml(mod2)*100, 2)# heteroscadasticity modelmod2b <-rma.mv(yi = SMD, V = VCV, mod =~ Type -1, random =list(~1|FocalSpL , ~1| RecNo, ~Type | Obs_ID), struct ="DIAG",test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod2b)# make other heteromod2e <-rma.mv(yi = SMD,V = VCV,mod =~ Type, random =list(~1|FocalSpL , ~1| RecNo, ~Type | Obs_ID), struct ="DIAG",test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod2e)mod2f <-rma.mv(yi = SMD,V = VCV,mod =~relevel(Type, ref ="LifeHistory"), random =list(~1|FocalSpL , ~1| RecNo, ~Type | Obs_ID), struct ="DIAG",test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod2f)# heteroscadasticity model better than the homoscedasticity model# note LifeHistory has lowest variation but this may be expected? # as it is less flexiable (e.g. the number of eggs?)anova(mod2, mod2b)orchard_plot(mod2b, mod ="Type",group ="RecNo",xlab ="Standardised mean differnece (SMD)",branch.size =3)```### trait categories```{r}# Category of responsesmod3 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Category -1, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod3)round(r2_ml(mod3)*100, 2)orchard_plot(mod3, mod ="Category",group ="RecNo", xlab ="Standardised mean differnece (SMD)",angle =45,branch.size =3)```### Predactor guild```{r}# Predactor guild# quite heterogeneous# TODO this could be in random effects - think abou thtis a bit latermod6 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ PredGuild -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod6)round(r2_ml(mod6)*100, 2)orchard_plot(mod6, mod ="PredGuild",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Setting```{r}# Settingmod7 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Setting -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod7)round(r2_ml(mod7)*100, 2)orchard_plot(mod7, mod ="Setting",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Season```{r}# Seasonmod8 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Season -1,test ="t",method ="REML",sparse =TRUE,data = dat)mod8b <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Season,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod8b)round(r2_ml(mod8)*100, 2)orchard_plot(mod8,mod ="Season",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Design```{r}# Designmod9 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Design -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod9)mod9b <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Design,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod9b)round(r2_ml(mod9)*100, 2)orchard_plot(mod9,mod ="Design",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Response period```{r}# Response periodmod10 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ ResponsePeriod -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod10)round(r2_ml(mod10)*100, 2)orchard_plot(mod10,mod ="ResponsePeriod",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Control type```{r}# control typemod11 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ ControlType -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod11)round(r2_ml(mod11)*100, 2)orchard_plot(mod11,mod ="ControlType",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Sex```{r}# sex# TODO - this could be interesting# what is in males and femalesmod12 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Sex -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod12)round(r2_ml(mod12)*100, 2)orchard_plot(mod12,mod ="Sex",group ="RecNo",xlab ="Standardised mean differnece (SMD)")# shoter data for just males and females# hetero but no sex effect heredat_sex <- dat %>%filter(Sex !="both")VCV3 <-vcalc(vi = dat_sex$Vd,cluster = dat_sex$SubjectID,rho =0.5)mod12a <-rma.mv(yi = SMD, V = VCV3, mod =~ Sex, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), #struct = "DIAG",test ="t",method ="REML", sparse =TRUE,data = dat_sex)mod12b <-rma.mv(yi = SMD, V = VCV3, mod =~ Sex, random =list(~1|FocalSpL , ~1| RecNo, ~Sex | Obs_ID), struct ="DIAG",test ="t",method ="REML", sparse =TRUE,data = dat_sex)summary(mod12b)anova(mod12a, mod12b)orchard_plot(mod12b,mod ="Sex",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Age```{r}# agemod13 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Age -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod13)round(r2_ml(mod13)*100, 2)orchard_plot(mod13,mod ="Age",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Predactor type```{r}# type of predatordat$PredTo[dat$PredTo ==""] <-NAmod14 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ PredTo -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod14)round(r2_ml(mod14)*100, 2)orchard_plot(mod14,mod ="PredTo",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Treatment duration```{r}# treatment durationdat$ln_duration <-log(dat$duration_days)mod15 <-rma.mv(yi = SMD,V = VCV,random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mods =~ ln_duration,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod15)mod16 <-rma.mv(yi = SMD,V = VCV,random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mods =~ ln_duration*Type,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod16)round(r2_ml(mod15)*100, 2)bubble_plot(mod15,mod ="ln_duration",group ="RecNo",xlab ="log(duration in days)",g =TRUE) +geom_point(data = dat,aes(x = ln_duration, y = SMD,color = Type,fill = Type,size =1/sqrt(Vd)), alpha =0.6) +scale_color_discrete() +#+ # how to put the legend for colourguides(color ="legend")#p + geom_point(aes(colour = Type))#scale_colour_manual(values = c("red", "blue", "green"))```## Meta-regression: multi-moderator (mostly)### full models```{r}######################## Mulit-variable models#######################dat_short$sln_duration <-scale(log(dat_short$duration_days))mod_full <-rma.mv(yi = SMD, V = VCVs, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~#Design + sln_duration*Type + sln_duration*Treat_mod + Sex,test ="t",method ="REML",sparse =TRUE,data = dat_short)summary(mod_full)round(r2_ml(mod_full)*100, 2)orchard_plot(mod_full,mod ="Type",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)orchard_plot(mod_full,mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)int_type <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Type")bubble_plot(int_type, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3)int_trt <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Treat_mod")bubble_plot(int_trt, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3)``````{r}# mulit-model selectioncandidates <-dredge(mod_full, trace =2)# displays delta AICc <2candidates_aic2 <-subset(candidates, delta <5) # model averagingmr_averaged_aic2 <-summary(model.avg(candidates, delta <5)) # relative importance of each predictor for all the modelsimportance <-sw(candidates)```## Publication bias### Funnel plot: uni-moderator```{r}```### Funnel plot: multi-moderator```{r}funnel(mod0r, yaxis="seinv",type ="rstudent")``````{r}funnel(mod_full, yaxis="seinv",type ="rstudent")```### Egger regression: uni-moderator```{r}# Eggerdat$effectN <- (dat$Ncontrol * dat$NTreat) / (dat$Ncontrol + dat$NTreat)dat$sqeffectN <-sqrt(dat$effectN)mod0e <-rma.mv(yi = SMD,V = VCV,mods =~ sqeffectN,random =list(#~1 | Phylo,~1| FocalSpL,~1| RecNo,#~1 | SubjectID, # incoprated as VCV~1| Obs_ID),#R = list(Phylo = cor_tree),test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod0e)bubble_plot(mod0e,mod ="sqeffectN",group ="RecNo",xlab ="Effective N",g =TRUE)```### Decline effect: uni-moderator```{r}# decline effectmod0d <-rma.mv(yi = SMD,V = VCV,mods =~ Year,random =list(#~1 | Phylo,~1| FocalSpL,~1| RecNo,#~1 | SubjectID, # incoprated as VCV~1| Obs_ID),#R = list(Phylo = cor_tree),test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod0d)bubble_plot(mod0d,mod ="Year",group ="RecNo",xlab ="Publication year",g =TRUE)```### All together```{r}# full modeldat_short$effectN <- (dat_short$Ncontrol * dat_short$NTreat) / (dat_short$Ncontrol + dat_short$NTreat)dat_short$sqeffectN <-sqrt(dat_short$effectN)mod_fulle <-rma.mv(yi = SMD, V = VCVs, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ sqeffectN + Year + sln_duration*Type + sln_duration*Treat_mod + Sex,test ="t",method ="REML",sparse =TRUE,data = dat_short)summary(mod_fulle)dat_fulle <-qdrg(object = mod_fulle, data = dat_short)# marginalized overall mean at vi = 0 and year.c = 0; also weights = "prop" or "cells" average things over proportionally. if not specified, all groups (levels) get the same weights# res_fulle1 <- emmeans(dat_fulle, # specs = ~ sqeffectN,# df = mod_fulle$ddf, # weights = "prop")```